Work fluctuations for a Brownian particle in a harmonic trap with fluctuating locations 
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We consider a Brownian particle in a harmonic trap. The location of the trap is modulated according to an 
Ornstein-Uhlenbeck process. We investigate the fluctuation of the work done by the modulated trap on the 
Brownian particle in a given time interval in the steady state. We compute the large deviation as well as the 
complete asymptotic form of the probability density function of the work done. The theoretical asymptotic 
forms of the probability density function are in very good agreement with the numerics. We also discuss the 
validity of the fluctuation theorem for this system. 
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I. INTRODUCTION 

Equilibrium statistical mechanics provides us a well- 
established framework to deal with systems in thermal equi- 
librium. When a system is perturbed externally from its equi- 
librium state, the so-called fluctuation-dissipation theorem re- 
lates the linear response of the system (to the external per- 
turbation) to the fluctuations properties of the system in equi- 
librium (in the absence of the perturbation) [1|. Within the 
framework of linear response theory, the Green-Kubo rela- 
tion gives linear transport coefficients in terms of integral over 
time-correlation function of the corresponding current in equi- 
librium [2J. In contrast, a general understanding of nonequi- 
librium (arbitrarily far from equilibrium) systems is rather 
poor. That is why there has been a lot of excitement surround- 
ing the fluctuation theorem, which aims at making a general 
statement about the fluctuations of entropy production during 
a nonequilibrium process. The fluctuation theorem has been 
suggested as a natural extension of the fluctuation-dissipation 
theorem from the linear response regime to arbitrarily far from 
equilibrium, as the fluctuation theorem reduces to the Green- 
Kubo formula and the Onsager reciprocity relations in the zero 
forcing limit . 

Imagine, two heat reservoirs at different temperatures con- 
nected by a thermal conductor. For a macroscopic object, we 
expect the heat to flow across the conductor from the hot- 
ter to the colder reservoir, in accordance with the second law 
of thermodynamics. However, for small systems, where mi- 
croscopic fluctuations become important, once in a while we 
might observe heat to flow from the colder to the hotter reser- 
voir. These reverse events are usually referred to as "second 
law violation". The fluctuation theorem gives a mathemati- 
cal expression for the ratio of the probability of "obeying the 
second law" to that of the "second law violation". 

Following a theoretical argument, Evans et. al. found 
a relation between the probabilities of positive and negative 
entropy production in the nonequilibrium steady state, in a 
molecular dynamics simulation of a two-dimensional fluid 
driven by external shear and coupled to a thermostat |4j. 
Gallavotti and Cohen proved this relation (and called it fluctu- 
ation theorem) for the phase space contraction (interpreted as 
the entropy production) for dissipative dynamical systems in 
the nonequilibrium steady-state using chaotic hypothesis and 
time reversal invariance Q. Evans and Searles had derived 



earlier a similar relation (now known as the transient fluctu- 
ation theorem) for systems starting from equilibrium initial 
condition (6). For stochastic systems, the fluctuation theo- 
rem has been proven by Kurchan Q for Langevin dynamics 
and extended by Lebowitz and Spohn (8) to general Markov 
processes. Subsequently, there has been an explosion of re- 
search activities investigating the validity of the fluctuation 
theorem for other quantities such as work, power flux, heat 
flow, total entropy, etc., both theoretically Il9j j26l and experi- 
mentally E71 - I4T1 . The recent review ll42l contains an exten- 
sive list of references pointing to several other reviews as well 
as research articles on fluctuation theorem and related topics. 

Recently, Ref. [39 1 reported experiments on the fluctuations 
of the work done by an external Gaussian random force on 
two different stochastic systems coupled to a thermal bath: (i) 
a colloidal particle in an optical trap and (ii) an atomic-force 
microscopy cantilever. Analytical results have been obtained 
for the second system in [24|. In the first experiment, a col- 
loidal particle immersed in water (which acts as thermal bath) 
is confined in an optical trap. The position of the trap is mod- 
ulated according to a Gaussian Ornstein-Uhlenbeck process. 
The authors have experimentally determined the probability 
density function (PDF) of the work done on the colloidal par- 
ticle by the random force exerted by the modulating trap. In 
this paper, we analytically treat this problem. 

The remainder of the paper is organized as follows. In 
Sec. |IlJ we define the model. Section |Hll contains the deriva- 
tion of the moment generating function of work done W T in a 
given time T, which has the form (e~^ Wt ^ g(X)e r ^^ for 
large T. In Sec. IV we invert the moment generating func- 
tion to obtain the asymptotic form (for large t) of the PDF of 
W T . We find that in that relevant interval, g(X) can either be 
analytic, or can have either one branch point or three or four 
branch points, depending on the values of the tuning parame- 
ters of the problem. The case when g (X) is analytic, is simpler 
and the asymptotic PDF can be obtained by the usual saddle 
point approximation, which is given by Eq. (|27]i in Sec. IV A 



In Sec. IV B we deal with case when g{X) has one branch 



point. The cases when g{X) has three and four branch points 
are discussed in Sees. |IVC| and IV D respectively. The an- 
alytical results obtained in each section are supported by nu- 
merical simulation performed on the system. Section [V] con- 
tains a discussion on large deviation function and validity of 
the fluctuation theorem in the context of the problem at hand. 
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Finally, we summarize the paper in Sec. VI Some details of 
calculation are pushed to two appendices: Appendix [A] con- 
tains the details of the calculation of the moment generating 
function. In Appendix |Bj we analyze the singularities of g(X) 
and in Appendix |C| we give the steepest descent method with 
branch point to calculate the PDF. We also provide an index 
of the relevant notations used in the paper in Appendix [P] for 
the ease of quick look up. 



II. THE MODEL 

Consider a Brownian particle suspended in a fluid at tem- 
perature T, with the viscous drag coefficient y. The particle 
is confined in a quadratic potential (harmonic trap) around the 
position y and having a stiffness k. The position x(t) of the 
particle is described by the overdamped Langevin equation 



dx 
dt 



y 



(i) 



where = y/k is the relaxation time of the harmonic trap. 
The thermal noise 4(f) is taken to be Gaussian with mean 
(4(f)) = and covariance (£(t)t;(s)) = 2DS(t — s), where 
the diffusion coefficient D = y~ l ksT with kg being the Boltz- 
mann constant. An external time-varying random force is 
exerted by the trap on the Brownian particle by externally 
modulating the position of the trap according to an Ornstein- 
Uhlenbeck process 



dy 
dt 



y_ 

TO 



C(0 



(2) 



where £(f) is an externally generated Gaussian white 
(non-thermal) noise with mean (4(f)) = and covariance 
(£(f )C( S )) =2AS(t — s). There is no correlation between the 
externally applied noise and the thermal noise, (C,(t)^{s)) — 
0. The system eventually reaches steady state, and in the 
steady state the trap exerts a correlated random force ky(t) 
on the Brownian particle with mean (>'(f )) = and covariance 
(y(t)y(s)) — AToexp(— \t — s\/tq). The quantity of our inter- 
est is the work done in the steady state, by the random force 
exerted by the trap on the Brownian particle in a given time 
duration T. This is given (in units of kgT) by 



1 f T dx 
W T = — / ky{t)—dt, 
kgl Jo dt 



(3) 



with the initial condition (at 1 = 0) drawn from the steady 
state distribution. 

It is convenient to use the following dimensionless param- 
eters 



0=A/D, and 8 = Tq/t- 



(4) 



From an experimental perspective [39], it is natural to use 
another parameter that measures the deviation of the system 
from equilibrium: 



where (x 2 ) is the variance of x in the steady state in the pres- 
ence of trap modulation, whereas (x 2 ) eq = D%y is the corre- 
sponding variance at equilibrium, i.e., without the presence of 
the trap modulation (y = 0). It should be noted that, the three 
parameters introduced above are not independent of each oth- 
ers and are related by 



a = e8 2 (l + 8) 



-1 



(6) 



The mean work can be computed easily using the above 
equations and one finds (W T ) « o:t/to for large T. Although 
the mean work is positive (and large for large t), there can be 
negative fluctuations (with small probabilities) and the fluc- 
tuation theorem quantifies the ratio of the probabilities of the 
positive and the negative fluctuations. However, the aim of 
this paper is not to merely check whether the fluctuation the- 
orem is satisfied or not for this system, but to obtain the dis- 
tribution of the work done (for large t), which contains much 
more information about the system than the former. 



III. MOMENT GENERATING FUNCTION 

To compute the distribution of W T , we first consider the mo- 
ment generating function restricting to fixed initial and final 
configurations (xo,yo) and (x,y) respectively: 

Z(X, x ,y,T\xo,yo) = (e- Xw *8[x-x(T)]8\y-y(z)]) {xoyo) , 

(7) 

where (• • • )(x ,yo) denotes an average over the histories of 
the thermal noises starting from the initial condition (xo,yo). 
It can be shown that Z(X,x,y,T\xo,yo) satisfies the Fokker- 
Planck equation 



dz 

dx 



(8) 



with the initial condition Z(X,x,y,Q\xo,yo) — 8(x~ xq) 8(y- 
yo), and the Fokker-Planck operator is given by 



dx 2 



„ d 2 8 d , . Id 
ay L To dx To ay 

X8 2 , . X 2 8 2 



2X8 d 

Tq " dx 



+ — y ^ + T 2 B yiX - y) + ^ r 



(9) 



We do not know whether the above partial differential equa- 
tion can be solved to obtain Z. Fortunately, however, one does 
not require the complete solution of the above equation to de- 
termine the large-T behavior of the distribution of W T . 

The solution of the Fokker-Planck equation can be formally 
expressed in the eigenbases of the operator and the large 
T behavior is dominated by the term having the largest eigen- 
value. Thus, for large-T, 



Z{X,x,y i x\ X Q,yo) = x(x ,y ,X) x i J (x,y,X)e 



, (10) 



(5) 



where *P(x,y,X) is the eigenfunction corresponding to the 
largest eigenvalue fl(X) and %{xo,yo,X) is the projection of 
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the initial state onto the eigenstate corresponding to the eigen- 
value jtt(A). While we cannot solve the Fokker- Planck equa- 
tion, the functions in Eq. ( [T0| can be obtained using a method 
developed [22 1 and used |24] recently. We sketch the deriva- 
tion in the context of the present problem in Appendix [A] 
where we find 



M(A) = ~[l-v(A)] 
2 re- 
in which v(A) is given by, 



T c = T (l + S)-\ (11) 



v(A) = y/l+4ak(l-k), a = a(l + Sy 



(12) 



We observe that the eigenvalue satisfies the Gallavotti-Cohen 
symmetry /i(A)=/i(l— A). In terms of the column vector 
U = (x,y) T , the eigenfunctions are 



«P(jt,y,A) 



1 



X{x 0l y ,k) =exp 



2a;^/det//i(A) 
1 



exp 



- l -u T u{X)u 



UiL 2 {X)U Q 



(13) 
(14) 



where the matrices H\,L\, and L2 are given in Appendix [A] 

Using the explicit forms one can verify the eigen- 
value equation Jzf^^Xjy, A) = /J.(X) l ¥(x,y, A). Moreover, 
J-ooJ-ooXixiyi^^iXjyj^dxdy — 1, as expected. From 
the above expressions, we also find that jj.(Q) — and 
X(xo,yo,0) = 1. Since the A — case of Eq. (|7]i gives the PDF 
of the variables (x,y ) and ji (0) is the largest eigenvalue, it fol- 
lows from Eq. ( 10 1 that ^(xjjO) is the steady-state PDF of 



(x,y). Therefore, averaging over the initial variables (xo,yo) 
with respect to the steady-state PDF *¥(xo,yo>0) and integrat- 
ing over the final variables (x,y), we find the moment gener- 
ating function of the work in the steady state as 



Z(X,T) = {e- m <)=g(X)e^ + 



(15) 



where 
with 



v / v(A) + 1 - 2b+k a/v(A) + 1 - 2Z?_A 
2v(A) 

yjv{k) + l+2b+k yjv(k) + l+2b_X ' 



b± = 



a 



1±W1 



4 

65 



(16) 



(17) 



The first factor in the above expression of g(k) is due to the 
averaging over the initial conditions with respect to the steady- 
state distribution and the second factor is due to the integrating 
out of the final degrees of freedom. 



IV. PROBABILITY DENSITY FUNCTION 

The PDF of the work done W T can be obtained from the 
moment generating function Z(A,t), by taking the inverse 



"Fourier" (two-sided Laplace) transform 

1 r+i 



P(W r ) : 



2ni J-i 



Z(k,T)e lw *dk, 



(18) 



where the integration is done along the imaginary axis in the 
complex A plane. Using the large-T form of Z(A , t) given by 
Eq. ( 15 1 we write 



P(W T = wx) 



1 

2%i J-i 



+100 



g(k)e tfwW dk, 



where 



f w (k) = -[l-v(k)]+kw. 



(19) 



(20) 



and we have set T c = 1 for convenience. This is completely 
equivalent to measuring the time in the unit of z c , that is, 



The large-T form of P(W Z ) can be obtained from Eq. ( 19 1 
by using the method of steepest descent. The saddle-point /P 
is obtained from the solution of the condition f w {k*) = as 



k*(w) = 



Vw 2 +( 



(21) 



From the above expression one finds that k*(w) is a mono- 
tonically decreasing function of w and k*(w —> =F°o) — > k±, 
where 



A± = 



1 



1±W1 



(22) 



Therefore, A* G (A_,A + ). It is also useful to note that v(A) 
can be written in terms of k± as 



v(A) = v/4fl(A+-A)(A-A_). 



(23) 



This clearly shows that v(A) has two branch points on the 
real -A line at A±. However, v(A) is real and positive in the 
(real) interval A £ (A_ , A+ ) . As a consequence, f w (A ) remains 
real in the interval (A_,A+). At A = A* we find 



v(A*) = 



and 



h s (w):=f w {k*) = 
One also finds that 



1 + w — \J w 2 + a^\ 



2(w 2 +a) 3 / 2 
y/a(l+a) 



>0. 



(24) 



(25) 



(26) 



This means that f w (k) has a minimum at A* along real-A, 
and hence the path of steepest descent is perpendicular to the 
real-A axis at A = A*. 
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FIG. 1. (Color online) The regions in the (a) 6,8 and (b) a, S spaces, where g(X) has the number of singularities mentioned in the figure. The 
equations of the boundary lines separating different regions are given in Appendix|B] Ofc, =1/3 for 5 = and 6 Cn — > 1/3 as 5 — > °°. Each of 
the phase boundaries meet at 6 = (a = 0), 5 = 1. 



Now, if g(X) is analytic for X € (Q,X*), one can deform 
the contour along the path of the steepest descent through the 
saddle-point, and obtain P(W T ) using the usual saddle-point 
approximation method. However, if g(X) has any singulari- 
ties, then the straightforward saddle-point method cannot be 
used, and one would require more sophisticated methods to 
obtain the asymptotic form of P(W T ). Therefore, it is essen- 
tial to analyze g(X) for possible singularities. In Appendix [B] 
we examine the terms under the four square roots in the de- 
nominator of g(X) in Eq. ( 16 1. 



In Fig. [T] we show the regions in the (8,5) and (a, 5) 
planes, where g(X) possesses singularities. 



A. The case of no singularities 




-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 



w 



In the singularity free region 5 < 1, 9 < 9 Cl (a < a Cl ), the 
asymptotic PDF of the work done is obtained by following the 
usual saddle-point approximation method. We get 



FIG. 2. (Color online). P(W r ) against the scaled variable w = W T /t 
for T = 10, T c = 1. The points (blue) are obtained from numerical 
simulation, and the dashed solid lines (red) plot the analytical asymp- 
totic forms given by Eq. {27}. 6 Cl = 9/35 = 0.257 ... for 5 = 1/2. 



p{yv x = wx) 



(27) 



where h s (w) and /"(A*) are given by Eqs. ( |25] > an d p6| ), re 



spectively, and g(X*) can be obtained from Eq. (I16J while 
using X* from Eq. \2\\ . Figure |2| sh ows very good agreement 
between the form given by Eq. (27 i and numerical simulation 
results for 9 < 9 C{ . 



B. The case of one singularity 

In the case 8 < 1, 9 > 9 C] where g(X) has only one singu- 
larity, or the case 5 = 1 where only one singularity of g(X) is 
relevant, we can write 



gi(A) ^ 

y/X a — X 



(28) 
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where g\ (A) is the analytical factor of g(X). 

It is evident that for a given value of 8 and 9, the position 
of the branch point A a is fixed somewhere between the origin 
and A+. On the other hand, according to Eq. (21 1, even for a 



where 



fixed 9, the saddle-point A*(w) moves unidirectionally along 
the real-A line from A_ to X + as one decreases w from +°° to 
— °° in a monotonic manner. Therefore, for sufficiently large 
w, the saddle-point lies in the interval (A_, A a ), and therefore, 
the contour of integration in Eq. ( |T9] > can be deformed into the 
steepest descent path (that passes through X*) without touch- 
ing A a (see Fig. [9J. However, as one decreases w, the saddle- 
point hits the branch-point, A*(w a ) = A a , at some specific 
value w — wl given by Eq. (Bill. For w < w a , since A* > A a , 



the steepest descent contour wraps around the branch-cut be- 
tween A a and A* as shown in Fig. 10 Leaving the details of 
the calculation to Appendix [C] here we present the main re- 
sults. 



1. w > w* 



For w > w*, following Appendix C 1 we get 

g (X*yh s (w) 



1 ; v^t/^A*) V ' V 

where the function R\ (z) is given by 

R l (z) = ^=^ 2 K 1/4 (z 2 /2), 



(29) 
(30) 



with Kiu(z) being the modified Bessel function of the sec- 
ond kind. It follows from the asymptotic form of Ki^(z) 
that Ri(z ->• °°) ->• 1. Therefore, for w > w*, Eq. |29] i 
approaches the form of the usual saddle-point approxi- 
mation given by Eq. ( [27) 1. On the other hand, using 
K 1/4 (z) ~ (l/2)r(l/4)(z/2)- 1 / 4 for small z, we get Ri (z) ~ 
r(l/4)y / z/2^. As w — > w* from above, i.e., when the sad- 
dle point approaches the branch point from below, h a (w) — 
h s (w) = f w {K)-f w {^) ^ (A a - A*) 2 /"(A*)/2. Therefore, 
the expression given by Eq. p9| ) remains finite, even when the 
saddle point approaches the singularity, i.e., 



P(W Z = wt) 



T(l/4) gl (X*)e 
2k [2t/»(A*)]V4 



as w — > vv a 



(31) 



«3(2) 



*2(Z), 



(34) 



°°. Therefore, 



with /?2(z) being given by Eq. (C20i. Using the asymptotic 
forms of Ri{z) given in Appendix C 2 a we get R^ (z) — >• 1 in 
the limit z 



y/x*\fM\ 



for w<w* 



(35) 



As w — > w* (from below), Pq (w, t) — > 0. 

The contribution coming from the saddle point is given by 
(see Appendix |C2b[ ), 



P s (w,x) 



,Ths(w) 



R 4 Ui;[h a (w)-h s {w)] 



where the function /?4 (z) is given by 

^(z)= 1 /fz/ /2 [/- 1/4 (z 2 /2)+/i/4(z 2 /2)] 



2j F 2 (1/2,1;3/4,5/4;z 2 



(36) 



4z 



(37) 



where 7 ±1 / 4 (z) are modified Bessel functions of the first 
kind and 2Fi{a\ ,ai',b\ 7 b2\z) is the generalized hypergeomet- 
ric function, defined by Eq. (C29 1. The small and large z be- 



haviors of Rn{z) are given in Appendix C 2 b 

For w « w* we get P§ (w, t) <C Pr (w, t). On the other hand 



Ps (w, t) acquires the same limiting form as in Eq. ( 3 1 1, when 
w — > w* (from below). 



3. Numerical Simulation 

We now compare the asymptotic forms presented in this 
subsection with numerical simulation. In one case, we choose 
8 = 1 and 9 = 4, for which we get A± = (1 ± y/2) /2, A* (w) = 
(1 - v^w/Vl+w 2 ) /2, A a = 1/2, and w* = 0. In an an- 
other case, we choose 5 = 1/2 and 9 = 13.5, for which 
wl = —0.0135 Figure [3] shows very good agreement be- 
tween the analytical and and simulation results. 



C. The case of three singularities 



w <w„ 



For w < w*, following Appendix |C2| we write 

P(W z = wt) wP b (w,t)+P s (w,t), (32) 

where Pq (w, t) is the contribution coming from the integra- 
tions along the branch cut and P^(w,t) is the saddle point 



contribution. Following Appendix C 2 a we get, 



Mw,x)^^^=R^[^x[K{w)-h s {w)\ ). .33) 



Now we consider the case, 8 > 1 and 9 > 9 Co , in which 
case g(X) has three singularities (see Fig. [TJ at A a , Ac and 
Ad given by Eqs. jBl\ , ( TB3[ ) and ( |B4| i respectively; where 
A_ < Ac < < Aa < Ad < A+. Therefore, g(X) can be written 
as 



\/A — AcVAa — X\J 'Ad — A 



(38) 



where #3 (A) is the analytical factor of g(A). We notice from 
Eq. ( pTj i that A* — > A_ as w — > +°° and A* increases mono- 
tonically towards A + with decreasing w. Therefore, there are 
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FIG. 3. (Color online). P(W r ) against the scaled variable w = W x /x 
for T = 10, T c = 1. The points (blue) are obtained from numeri- 
cal simulation, and the dashed solid lines (red) plot the analytical 
asymptotic forms given by Eq. ( |29| for w > w* and Eqs. <[32[> — (|3<5[> 
for w < w a , where w a = for 5 = 1, 0=4, and w a = —0.0135 . . . 
for 5 = 0.5,6 = 13.5. 




FIG. 4. (Color online) Schematic steepest descent contours for the 
case when there are three branch points at A a , Ac and where A_ < 
Ac < < A a < A(j < A+; and the saddle point A* lies between (a) A_ 
and A c , (b) A c and Aa, (c) Aa and A^, and (d) A^ and A+ respectively. 



specific values +°° > w* > w* > w* d > —°° of w given by 
Eq. ( |B11 1 at which the saddle point hits the corresponding 
branch point, i.e., A*(w*) = A c , A*(w a ) = A a and X*(w* d ) = 
Ad- 



1. w> w* 



For w > w*, the saddle point lies between A_ and A c . Therefore, as in the case of one singularity discussed above in Sec. IV B 
the contributions comes from the branch point as well as from the saddle point, as shown in Fig.|4](a). Following the procedure 
similar to that in the one singularity case (see Appendix|C 2), we get 



P(W T = wt) 



y/™\fM*)\ 
|g(A*)k^W 



zh c (w) / , \ 



R 5 ( z[h c (w) - h s {w)] , ^T[h a (w)-h s (w)], \J*[h d (w)-h s (w)] 



where R3 (z) is given by Eq. ( 34 1, and 



r. I ZIZ2Z3 / _„2 

«5(Zl)Z2,Z3j = \l — — J due 



\Jz\ + iu\Jzi + iuy/z3 + iu \fz\ — iuy/zi — iu^/zs- 



(39) 



(40) 



2. w* < w < wJ 



For w* < w < w*, the saddle point lies between A c and A a . Therefore, the contour of integration can be deformed through the 
saddle point without crossing any singularity, as shown in Fig. [4] (b). Now, to compute the saddle point contribution one can 
follow the methods of Appendix |C 1| while taking into account of both the singularities A a and Ac. The calculation yields 



P{W T = wt) 



g(k*)e z H >( w 



==/? 6 ^t[h c (w)-h s (w)], ^T[A a (w)-A s (w)], yjz[h d {w)-h s (w)] 



(41) 



where 



^6(Z1,Z2,Z3) 



Z1Z2Z3 



e " du 



Vzi + iu\Jzi — iu^/z3 — iu 



(42) 
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As w — > w*, the first term of Eq. ( 39 1, coming from the integral along the branch cut, goes to zero. On the other hand, it can be 
shown that R${z\ — > 0, 22,23) = Rf>\Z\ — > 0, 22,23)- Therefore, Eqs. <j39j and (|4T]i approach the same limiting form as w —> w* c 
from the two sides. 



3. w* d < w < w* 



For w* A < w < w*, the saddle point lies between Aa and Ad- Therefore, the deformed contour is as shown in Fig. |4](c). 
Combining the contributions from the branch point A a and the saddle point, we get 



where 



P{W T = wt) 



,TAa(w) 



Ri 



Aa(w) - h s (w)] , 1/ T [Ad(w) - /z s (w)] 



loCl^le^W / / / / 

Xjxf"{X*) R% [y T M W )- h s( W )} > V T[*aH - As W] , ^[AdW-A.W] 



^7(21,22) 



2zi(z\+Z2) f Zl e 



,—2ziU+u 



: du, 



UWZl +Z2 — U 



and R s l;.\.:2.:\> = \/ due 1,2 



\JZ\ + M\JZ2 + iU\Jz-i — iu \pZ\ — IU\JZ2 — iuy/zi + » 



(43) 

(44) 
(45) 



As w — > w*, the first term of Eq. (43 1, coming from the integral along the branch cut, goes to zero. On the other hand, it can be 
shown that Re{z\,Z2 — > 0,23) = RgJii,Z2 — > 0,23). Therefore, Eqs. ( |4"T| i and (j43j approach the same limiting form as w —> w* 
from the two sides. 



W < Wj 



Finally, for w < w* d , the saddle point lies between Ad and A+. In this case, the integral along the branch cut can be divided 
into two parts: one, from Aa to Ad and another from Ad to A*. Between Ad and A*, the the integral above the branch cut exactly 
cancels the integral below the branch cut. Therefore, the net contribution is the sum of the contributions coming from the integral 
around the branch cut between A a and Ad, and the contribution of the integral along the contour (C\ and C4) through the saddle 



point, for which the calculation is similar to the one given in Appendix C 1 Therefore, we get 



where 



eTA \e xh '^ w > ( I / 

p(W r = wt) a R 9 J%[K{w) - h s (w)] , Jz[h d (w) - h s (w)] 



T/! S (IV) 



Rw[ \/t[hc(w)-h s (w)], Jr[h a (w)-h s {w)],Jt[h d {w)-h s (w)] , 



(46) 



#9(21,22) = 



2zi(zi-z 2 ) f Z1 ^ e- 2 ^ u+u2 du 



\/Uy/zi ~ Z2 ■ 



(47) 



#10(21, 22, 23) = 



212223 



% 



du 



-00 yjz\ + lUy / Z2 + iUy/Z3 + IU 

(48) 



It is evident from the above equations that Rj(zi,0) = 
Rg(zi,0). Moreover, it can be shown that #10(21,22,23 — >■ 0) = 
—R%(zi ,Z2,Z3 —> 0). Therefore, Eqs. ( |43| ) and |46) approach 
the same limiting form astv^yvj from the two sides. 



D. The case of four singularities 

Finally, we consider the case 8 > 1 and 9 < C ■, , in which 
case g(A) has four singularities (see Fig. [I} at A a , At,, A c and 
Ad given by Eqs. ( |B1| i-( |B4] i respectively; where A_ < Ab < 
A c < < A a < Ad < A + . Therefore, and g(X) can be written 

as 



8a{X) 



1/A — AbVA — A c %/ A a — A 



(49) 



where #4 (A) is the analytical factor of g{X). 

Now as w varies from +°° to — °°, the saddle point hits the 
branch points, A*(w*) = A; with i E {b, c, a, d}, at specific 



8 



values of w given by Eq. (Bill and 



-oo > vvjj > w r 



> w! > 



> — oo. It is straightforward to generalize the above results 



to this case of four singularities. Therefore, we only give the 
results below, without repeating the details. 



1. w > wt 



For w > wj, the saddle point lies between A_ and A|,, and 



P(W Z = wt) 



Rg ( ^t[A c (w)-A 8 (w)], yT[/l b (w)-/l s (w)] 



ig^^^ ^n [A b (w) - A s (w)] , ^ T [A c (w) - A s (w)] , yj X [A a (w) - h s (w)] , yjx[h d (w) ~ h s (w)] ) , (50) 



where Rg{z\,Z2) is given by Eq. (47 i and 



flll(Zl,Z2, Z3,Z4) 



21222324 



e " du 



K J- 00 \Jz\ - lUy/z2 — lUy/zi — iuy/ ZA ~ I 



21222324 



e " du 



7t J-oc y/zi+ iu^/zi + iu^/z-i + iu^fzA + iu 

(51) 



2. wt < w < wt 



For w* < w < w£, the saddle point lies between Ab and A c , and 



f>(W T = wt) 



Tftc(vc) 



#7 f yjx[h c (w)-h s (w)], ^ x[h b (w) - h s (w)] 



+ !^1^! /; 12 [yfx [h h {w) - hjw)] , ^ T [/l c (w) - A s (w)] , yj X [K{w) - A s (w)] , ^ x[h d ( W ) - hjw)] ) , (52) 



where Rj (21,22) is given by Eq. (44 1 and 

fZ]Z2ZsZ4 



^12(21,22,23,24) 



7C Jo 



due " 



1 



\fz\ — IU\/Z2 + IU\/Z3 + lUy/ZA + tU y/zi + iu^J Z2 — iu^/z3 — iu^J z\ — iu 



i. (53) 



3. w* < w < w* 

For w* < w < w*, the saddle point lies between Ac and A a , and the PDF is given by 

*\ th s {w) 



P(W r = wt) a ^2/"(aV 13 (V^bM-AsM], ^[AcH-AsH], v /t[A 3 (w)-A s (w)], ^vHW]) , (54) 



where 



„2 



P C \ /21222324 f e " c/m 

^13(21,22,23,24) = \ — - — / 1 — 1 — =f , —, (55) 

7t J -00 y/Zi +lUy/Z2 + lU^Z3 — lUy/Z4 — IU 



4. w* d <w <Wq 
For w* A < w < w*, the saddle point lies between A a and Ad, and 
P(^ T = wt) a Jg 7 Jt [A a (w) - A s (w) , Jt [A d (w) - A s (w) 

+ (\A ~ ' V T [ /tc(vv) ~ ' V T N w) - W] ■ \/ T [^(w) - As W] j , (56) 
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where R-j (21,22) is given by Eq. (44 1 and 



Ru{zi,Z2,Z3,Z4) = \/ ZlZ2 ^ Z4 J due' 



\Jz\ + iuy/z2 + iuy/zi + iuy/iA — iu y/zi — iu^/za. — iu^/z3 — Ju\/za + i 



i. (57) 



W < W j 



Finally, for w < w% the saddle point lies between Ad and A + , and 

«(A»)e**»M 



P(W T = WX) : 



/? 9 



/j a (w) -h s (w)],Jx [h d (w) - h s (w)] 



==Ru ( sj t Mw) - A s (w)] , a/ T [Ac(w) - h* (w)] , ^ X [A a ( w ) - A S ( W )] , yj x [h d (w) - A, W] J . (58) 



|g(A*)|e^W 

where #9(21,22) and #11(21,22,23,24) are given by Eqs. $7) and ( |5T] >, respectively. 

I 



It can be shown that, when w — > w* with i E {a,b,c,d} from 
the two sides of w*, the respective expressions of the PDFs, 
i.e, Eqs. (BO) and Eqs. (|52j and Q, Eqs. Q and Q , 
and Eqs. ( |56| l and ( |58| l, respectively, approach the same limit- 
ing form. 



6. Numerical simulation 

We now compare the analytical results obtained in this sec- 
tion with numerical simulation. We consider 5=5, for which 
we have 9 C2 = 0.18. Therefore, g(X) has three singularities 
for > C2 , whereas for 9 < C2 it has four singularities. 

For = 0.5 the three singularities are located at A c = 
-0.3062..., A a = 0.3958... and X d = 1.3062..., whereas 
A_ = -0.4848 ... and A+ = 1 .4848 .... Figure [5] compares 
numerical simulation for this case with analytical results ob- 
tained above for the case of the three singularities. 

On the other hand, for 6 =0.1, the four singularities of 
g(X) are located at at Ad = —10/7, A c = —1, A a = 13/11, and 
A d = 2. Moreover, A_ = -1.4621 ... and X+ = 2.4621 .... 
Figure |6]compares numerical simulation for this case with an- 
alytical results obtained above for the case of the four singu- 
larities. 
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FIG. 5. (Color online). P(W t ) against the scaled variable w = W t /x 
for T = 10, % c = 1, 5 = 5, and = 0.5. The points (blue) are obtained 
from numerical simulation, and the dashed solid line (red) plots the 
analytical asymptotic forms given in the text. The vertical dashed 
lines mark the positions w* = 0.8398 . . . , w a = 0.06269 . . . and wjj = 
-0.8398.... 



T . LARGE DEVIATION FUNCTION AND FLUCTUATION 
THEOREM 

The large deviation function is defined by 



h(w) = lim - lnP(W T = wx) 

T— >°o X 



(59) 



In other words, the large deviation form of the PDF refers to 
the ultimate asymptotic form P(W T = wx) ~ e T/i ( w ' while ig- 
noring the subleading corrections. Apart from being an inter- 
esting quantity on its own, the large deviation functions have 



found importance recently in the context of the fluctuation the- 
orem. The latter refers to the relation 



lim - In 

T— >oo X 



P(W T = +wx) 



P(W t 



-wx) 



w. 



(60) 



When the above relation is valid, the large deviation function 
evidently satisfies the symmetry relation 



h(w) — h{—w) = w. 



(61) 



Now, as we have seen in the above sections, when g(X) 
is analytic between the origin and the saddle point, the dom- 
inant contribution to P(W T ) comes from the saddle point as 
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the large deviation function is given by 



FIG. 6. (Color online). / > (W T ) against the scaled variable w = 
Wt/t for T = 10, T c = 1, 5 = 5 and = 0.1. The points (blue) are 
obtained from numerical simulation, and the dashed solid line (red) 
plots the analytical asymptotic forms given in the text. The vertical 



dashed lines mark the positions = 1 .4062 . 



= 0.3125, w* 



-0.0976 ... and w* A 



-0.3125. 



h(w) 



h c (w) for w > w* , 
h s (w) forw*<w<w* 
h a (w) forw<w*. 



(63) 



Since A c < 0, it is evident from Eq. (Bill that w* > 0. There- 
fore again, it is only when w* < (e.g., when 9 < 0.365 . . . 
for the 5 = 5 case), the symmetry relation Eq. (61 1 (and hence 
the fluctuation theorem) is satisfied only in the specific range 
max(w*,-w*) < w < min(-w*,w*). 

Therefore, for any 5, there exists a 9 C , given by w* = 
(equivalently A a — 1/2) as 



Q , s . 3 + 25 + 35 2 + (l-5)V9 + 145 + 95 2 

Be(S) = , (64) 

and the fluctuation theorem is not valid for 9 > 9 C . The 9 — 
9 C (8) line corresponds to the a = a c (8) line in the (a, 8) 
plane where 

, s . 3 + 25 + 35 2 + (l-5)\/9 + 145 + 95 2 
a c (8) = 2(TT ^ . (65) 

Figure [7] summarizes the state of validity of the fluctuation 
theorem in the 8, 9 and a, 8 parameter spaces. 



given by Eq. (27 1. On the other hand, when there are sin- 
gularities between the origin and the saddle point, the most 
dominant contribution to P(W T ) comes from the singularity 
closest to the origin (farthest from the saddle point) and lies 
between the origin and the saddle point. This is because, evi- 
dently — v(A), and hence the function f w (X), is convex on the 
interval [A_ , A+] and / W (A) is minimum at the saddle point A* 
along the real- A line. 

Consequently, for the case 8 < 1 and 9 < 8 Cl , where g(X) 
is analytic on the interval (A_,A+), the large deviation func- 
tion is h(w) = h s (w), given by Eq. (25 1. In this case, h(w) 
satisfies the above symmetry relation (|6l|>, and therefore, the 
fluctuation theorem is valid. On the other hand, for 8 < 1 
and 9 > 9 Cl , where g(X) has one singularity at A a , (also for 
8 = 1 and all values of 9, where only the singularity at A a is 
relevant), one has 



h(w) 



h s (w) for w > w* 
h a (w) for w < w* 



(62) 



Therefore, it is only when w* < (e.g., 
the 5 = 1 case), the symmetry relation Eq. (61 



when 9 < 4 for 
(and hence 



the fluctuation theorem) is satisfied only in the specific range 
w* < w < —w* A . Otherwise it is not satisfied. 

For the case 8 > 1, although there are either three or four 
singularities depending on whether 9 > 9 Cl or 9 < 9 C2 , the 
singularities closest to the origin (one on each side), namely A c 
and A a are common in both cases. Therefore, for both cases, 



VI. SUMMARY 

Let us now summarize the main contents of the paper. We 
have obtained analytical results for a system studied recently 
experimentally [39|. The experimental system consists of 
a colloidal particle in water and confined in an optical trap 
which is modulated according to an Ornstein-Uhlenbeck pro- 
cess. This system is described by a set of coupled Langevin 
equations. We have computed the PDF of the work done by 
the modulating trap on the Brownian particle in a given time 
T, for large z. The moment generating function of the work 
has the (e~^ Wr ) s» g(A)e T/J ^' for large T. Inverting this, we 
obtain the PDF of the work within the saddle point approxi- 
mation. 

The results can be described in terms of two independent 
parameters: (1) the parameter 9 that quantifies the relative 
strength of the external noise that generates the Ornstein- 
Uhlenbeck process for the trap modulation, with respect to the 
thermal fluctuations, and (2) the ratio 8 — Tq/Tj of the correla- 
tion time of the trap modulation to the viscous relaxation time 
of the particle in the trap without any modulation. We find that 
the cumulant generating function jU (A) is analytic in a (real) 
interval (A_,A+) and the saddle point lies within this interval 
— here A± depends on the values of 9,8. On the other hand, 
depending on the values of the pair (9,8), the function g(X) 
behaves differently. For 8 < 1, there exists a value 9 Cl (8) 
such that g(X) is analytic in the interval (A_, A+) for 9 < 9 C{ 
whereas it has a branch point for 9 > 9 Cl . For 8 > 1, there 
again exists a 9 C1 (8), and g(X) has either three or four branch 
points depending on whether 9 > 9 C2 or 9 < C , . For 5 = 1, 
there are three branch points of which two coincide with A± . 
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FIG. 7. (Color online). The phase diagrams in the parameters space showing the state of the validity of the fluctuation theorem. Apart from 
the 6 = 8 C and a = a c (orange) lines, the other lines are the same as in Fig.[TJ In the (light orange) regions above the (orange) lines 6 = G c in 
(a) and a = 0C C in (b), the fluctuation theorem is not valid at all, whereas it is always valid in the white regions below the (blue) lines 6 = 6 Cl 
in (a) and a = a Cl in (b). In the intermediate (light blue) region, the fluctuation theorem is valid only within a limited range of w given by 
w* <w < —w* for 8 < 1 and max(w*,-Wj) < w < min(— w^, w*) for 8 > 1. For 8 = 0, we have a Cl =1/3 and a c = 3. 



We have done the analysis in each of these regions and ob- 
tained the asymptotic form of the PDF accordingly. We have 
compared our analytical results with simulation results on this 
system and found very good agreement between the two. 

The calculation also gives the large deviation function as a 
by-product, using which we check the validity of the so-called 
fluctuation theorem for this context. We find that in the region 
8 < 1 , 9 < 9 Cl , it is always valid. Outside this parameter re- 
gion, there exists a 6 C (8) and the fluctuation theorem is valid 
for a limited range of w around zero when 9 < 9 C . For 9 > 9 C , 
the fluctuation theorem is not valid at all (see Fig.|7|. 

Finally, we would like to point out that it is very easy to 
generalize the analysis of the paper to the case of a Brown- 
ian particle subjected to any exponentially correlated external 
random force. Moreover, the asymptotic analysis (steepest 
descent method with branch points.) carried out in this paper 
should be applicable for finding asymptotic approximations of 
similar integrals in general. 



Appendix A: Calculation of the Moment Generating Function 

The evolution equations ([T| and |2) can be presented in the 
matrix form 



at To 



(Al) 



where U — (x,y) T and 7} = (E, , Q T are column vectors and A 
is a 2 x 2 matrix given by 



A 



8 8 
1 



(A2) 



Using the integral representation of 5-function, the re- 
stricted moment generating function defined by Eq. (|7]i 



Z(X,U,z\U ) = 



° T U / „-XW r -ia T U(x) 



{my 



(A3) 



where a T = (Gi,&2)- 

Substituting dx/dt form Eq. Q in Eq. ^ we get 



ACKNOWLEDGMENTS 



W T = / dt 
Jo 



-j-y(x-y) + —y% 



(A4) 
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which is useful to rewrite as 



1 f T 

W t = - dt 
I Jo 



^U T AiU + — {U T A T 2 T] +T] T A 2 U) 

T 



(A5) 
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where 



-1 
D V-l 2 



and A 2 = |(J JV (A6) 



Now, we proceed by defining the finite-time Fourier trans- 
forms and their inverses as follows: 

~ 1 f x 

[U(co„),ri(co n )} = - dt[U(t),rj(t)]exp(-iC0nt), (A7) 
X Jo 

P(f),T}(t)]= £ [^(o)„),^(co„)]exp(7c<v), (A8) 



with o„ = Innjx. In the frequency domain, the Gaussian 
noise configurations denoted by [77(f) : < t < t] is described 
by the infinite sequence \r\{(£>„) : n — — <»,...,— 1,0, +1, ...,°°] 
of Gaussian random variables with the correlation 



2D 



(ri(co)ri T (co')) = — 8(co + a>')diag(l, 0) 



(A9) 



In terms of the Fourier transform Eq. ( A5 1 can be written as 
1 



* = ? E 



2 -f/ r (a)„)Aif/(-co„) 



— {u T {(O n )A T 2 r]{-(o n ) + r] T {(o n )A 2 U{-(O n )\ 

To L J 



Equation ( Al 1 gives 



U((o n ) = i Gri(a> n )-^GAU, 

T 



(A10) 



(All) 



where G = [iuI + A]' 1 with u = (B„To, AU = £/(t) -t/(0), 
and / being the identity matrix. The elements of G are: 
Gn = (8 + iu)~ l , G22 = (1 + iu)~ 1 , Gu = 8GnG 2 2, and 
G21 = 0. Note that G(-u) = G*(u), 77 (-0)) = 7J*(a)), and 
{/(— £») = £/*(©). Substituting C7 in Eq. (A5l, and grouping 
the negative n indices together with their positive counterparts 
in the summation, we get 



T77 r (GjAi Go + A 2 G + G T qA t 2 ) 770 



E 

n=l 



-AfJ r (GjAiGo + G^)77o 
+ ^AC/ r (G r AiG )A{/ 

xi\ r {G T A l G* +A 2 G* +G T A 2 r )rj* 



~ ri T (G T A 1 G* +A 2 G*)AU 



-MJ T {G T A l G* +G T A T 2 )r]* 
1 

T 



At/ r (G r AiG*)Af/ 



(A12) 



in which Go = G(u = 0) = A 1 and 770 = 77 (0). 

Next, we we express U(x) in terms of the Fourier series 



E/(T) = lim V U{co„)e- ilo " e 

e->0 „_ „ 



(A13) 



While inserting U from Eq. ( Al 1 1 into the above equation, we 
observe that (1/f) Ge~' m " E — > as T — > 00. This is because 
in the large-T limit, the summation can be converted to an in- 
tegral which can be closed via the lower half plane, and the G 



is analytic there. Thus, using only the first term of Eq. ( Al 1 
we get 

o t U(t) = To cj 2 " Go 770 

+ t £ [e- ia, " £ rj T G T a + e i<0 " e o T G*rj*]. (A14) 

n=\ 

Using this expression as well as W T from above in Eq. ( |A3| >, 
we get 



Z(X,U,z\U ) = 



Jo' U 



where s„ is quadratic in 77, given by 



n<<*>. 



(A15) 



so 
and 



— ^Bom + aom-^^GlAiGoAU, (A16) 
2 IT 



s„ = - Xtt] 1 B„rj* + r} T a„ + a T n r]* 
X 

- -AU T G T AiG*AU for n>\, 

T 

in which we have used the following definitions: 

B„ = G T A\G* +A 2 G*+G T A T 2l 
Bo = G%A l Go+A 2 Go + G%Al, 
a„ = A [G r AiG* +A 2 G*]AU - i^e'^G 1 a, 
a T „ = XAU T [G r AiG* + G T A T 2 ] - iT e ia " £ a T G* . 



(A17) 



(A18) 
(A19) 
(A20) 
(A21) 



Therefore, calculating the average (e s ") independently for 
each n > 1 with respect to the Gaussian PDF P(r\) = 
7r" 2 (detA)- 1 exp(-77 r A- 1 77*) with A= (2D/r)diag(l, 0) 
we get 



(*>) 



exp(a T „Q.- l a n - ^AU T G T A\G*AU) 
det(A£2„) 



(A22) 



where £2„ = i(\B n + %~ 1 A - 1 ) . For the n = term, calculating 
the average (e s °) with respect to the Gaussian PDF P(r)o) — 
(2^) _1 (detA) _1 / 2 exp(-57j r A _1 7jo) we get 



^ = exp(ia r n V- ±- x AU T GlA x GoAU) 
K ' ^det(Ai2 ) 



Since, (e* n ) = (e s - n ), the product in Eq. ( A15 1 yields 



n=0 



n<^")=«P(-9 E ln[det(An„)] 



1 



11 = 



xexp — £ [a r „T£2- 1 a„-AAf7 7 'G 7 'AiG*Af/] 

, 2T „=_„, 



(A24) 



13 



The determinant in the above expression is found to be 



det(AQ„) = [l+40A(l-A)5 z M z |Gn| z |G 2 2r]. (A25) 

Now, taking the large z limit, we replace the summations over 
n by integrals over CO, i.e., — > tJj^. After, evaluating 
the integral, the argument of the exponential in first line of 
Eq. < |A2t| yields 



t/Tq 
47Z 



du In [det(A£2„)] = Tfi(X), 



(A26) 



where ju(A) is given by Eq. ( pT| . Similarly, converting the 
argument of the exponential in the second line of Eq. ( A24\ in 
the integral forms, after some manipulation we get 



: a T H 1 a+iAU T H 2 a- 



1 



; AW H 3 AU 
(A27)' 



in which Hi, Hi, and i/3 are given by 

£>T 



H\ = 



2n J- 



duG*Qr l G T 



H 2 = — lim — , 

e->0 2% J-c 



(A28) 

rfw e"' £ / T ° (G + A 1 G + G + A[ ) (fl- 1 ) * G f , 



(A29) 



A 2 1 r T~ T~T ~ 1 

H 3 = / du{G T A x G*+G T Al)Qr l 



x (G T AiG* +A 2 G*)- 



X 1 
27tDT 



[G r A[G*], 
(A30) 



where we have used where £2„ = z~ DQ. n and A] 2 = DAi 2 
so that the integrands remain dimensionless and dimensions 
are carried outside to the integrals. We then evaluate the in- 
tegrals performing the method of contours in the complex u 
plane, and using G* (it) + G(u) = 2GAG* and G* (a) - G(h) = 
2iuGG* , which yields 



#l(A): 

H 2 {X) 



Dtp 
8(l + 8)v(X) 

v(A)-l A 
2v(A) ^0 1 

Ag 2 
Dt (1 + <5)v(A) 



5 + 65 2 6<5 2 
05 2 98 + 08 

X8 fed 98 
1 



2 1 , (A31) 



(l + 5)v(A) 



X98 1 
1 A - 



(A32) 
(A33) 



Finally, inserting Eq. ( |A27| > in Eq. ( A15 1, and performing 
the Gaussian integral over a while using the facts that H\ and 
Ht, are symmetric and Hj, — H^ l H^ +// 2 //p 1 //J we get 



Z(X,U,t\U ) we T "W«p ( L 2 (X)U 



2^v / det7^(A) 



with Li(A) = fff 1 +H{ l H 2 T and L 2 (A) 



exp --I7 r Li(A)I7 , (A34) 



-H7 x Hl. From 



the above equation, it is trivial to identify %(Uo,X) and 



*F(t/,A) used in Eq. (IOi. Since, L\ +L 2 = H7 l , it is evident 
that fx{U,X)¥(U,XjdU = 1. 

Application of the Langevin operator given by Eq. Q on 
*P(t/,A) yields 



D(L 



1,1\2 



aD^' 2 ) 2 -^ 1 



x 2 y/(J/,A) 



+ 



1,2\2 



■CCDIL 



2,2 )2 . 



-(1-2A)L|' 2 



1 ,2.2 _ 
TO 



2 A(1-A) 



+ 



2DL}' 1 l[' 2 - 



-2DaL{' 2 L 2 ' 2 



-(1-2A)L 



5 rl,2 



Dt 2 ' 



—DL\ — aDL x 



y 2 v(u,X) 
8^ 

To 

xyy{U,X) 
5" 



T() 



VA(t/,A), (A35) 



where L^ 1 denotes the (/,y')-th element of the matrix L\. Us- 
ing the explicit expressions on the right-hand side of the 
above equation, after simplification, we find the coefficients 
of x 2 W(U,X), y 2x i>{U,X), and xyW(U,X) to be zero. The 
last term in square brackets in front of *P(E/,A) yields ju(A) 
given by Eq. (jTTJi. This verifies the eigenvalue equation 

j? x y(u,x) = [ni)y(u,x). 

The steady-state of the system is given by 



P SS {U)=V(U,0) 



expl 



U T H7 l {<d)U) 



27T\/det//i(0) 



(A36) 



Integrating Eq. \A3) over U and then averaging over the ini- 
tial condition Uq with respect to the steady state distribution 
Pss(Uq), we obtain Z(A) given by Eq. ( fl5| ), with 

g{X) = (det//i(A)det//i(0)detLi(A)det[//f 1 (0)+L 2 (A)])" 
= (det[l-v(A)i/ 2 r (A)]det[l+// 2 r (A)])" 1/2 , (A37) 

where to obtain the second expression, we have substituted 
the expressions of hi, L 2 and H\(Q) — v(X)H\(X). Inserting 
the matrix H 2 and evaluating the determinants, after simplifi- 



cation, we obtain Eq. ( 16 1 



Appendix B: Singularities of g(X) 



From Eq. (23 1 we recall that v(A±) = and v(A) > (is 
a semicircle) for A G (A_,A+). Moreover, all the four func- 
tions 1 ± 2b±X are linear in A with slopes ±2b± (where all 
four combinations of the two ± signs are considered). There- 
fore, for example, if (1 — 2b + X) has opposite signs at the two 
end points X±, then the function [v(A) + (1 — 2b + X)] must 
cross zero at some intermediate A. This is also true for the 
other three cases. From Eqs. ( fTTj i and |22} respectively, we 
note that b + > 0, b- < and A + > 0, A_ < 0. One can there- 
fore determine whether g (A) has a singularity as follows (see 
Fig.®: 
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(a) Evidently, 1 - 2b + X- > 0. Thus, v(A a ) + 1 - 2b+X a = 
for a specific A a G (A_ , A + ) if and only if 1 — 2b + X + < 0. 
When this happens [see Fig. [8] (a)], the position of the 
singularity can be found as 



K = {a + b+)/{a + b z + ). 
It is evident that Aa > 0. 



(Bl) 



(b) Evidently, 1 -2Z?_A+ > 0. Thus, v(A b ) + 1 -2Z?_A b = 
for a specific Ab € (A_ , A+) if and only if 1 — 2b- A_ < 0. 
When this happens [see Fig. [8] (b)], the position of the 
singularity can be found as 



X h = {a+b-)/{a + b±) 
and it can be shown that A b < 0. 



(B2) 



(c) Evidently, 1 + 2b+X+ > 0. Thus, v(A c ) + 1 + 2b+X c = 
for a specific A c € (A_, A + ) if and only if 1 + 2£> + A_ < 0. 
When this happens [see Fig. [8] (c)], the position of the 
singularity can be found as 



X c = {a-b+)/{a+b%) 
and it can be shown that A c < 0. 



(B3) 



(d) Evidently, 1 + 2£_A_ > 0. Thus, v(A d ) + 1 + 2b-X A = 
for a specific Ad € (A_ , A + ) if and only if 1 + 2b_A + < 0. 
When this happens [see Fig. [8] (d)], the position of the 
singularity can be found as 



X d = (a-b_)/(a + b 2 _). 



(B4) 



It is evident that Ad > 0. Moreover, it can be shown that 
A c + Ad = 1 . 

It is easily seen that the singularities of g(X) are branch points 
(square root singularities) and the function / H (A) at these sin- 
gularities is given by 

hi(w) := / V ,(A) = i [1 - v(A)] + A,-w, (B5) 

where the index i stands for one of the indices from the set 
{a, b, c, d}. Substituting v(A,) at the singularities using the 
conditions from above, we get 



h a (w 
K{w 
h c (w 



= (1 -fr + A a ) +A a w, 
= (1 -b-X b ) + A b w, 
= (1 +b + X c ) +A c w, 
= (1 +b-Xd) + XdW. 



(B6) 
(B7) 
(B8) 
(B9) 



It is also useful to define the non-singular part of g(X) at a 
singularity as 



Si^i) = lim |(A - A,y/ 2 g(A)| 



(B10) 



We note that the for a given set of parameters 9 (or a) and 
5, the position of the singularities (whenever they exist) are 
fixed within the interval (A_, A+). The specific values of w at 
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FIG. 8. (Color online) In the shaded regions of the (0,8) plane in the 
figures (a), (b), (c) and (d), the respective mathematical conditions 
given there are satisfied and consequently possesses singulari- 
ties at A a , Ab, Ac, and Ad respectively, given by Eqs. JbTI — ([B4ll. 



which the saddle point coincides with one of the singularities 
is obtained by solving A*(w*) = A; as 



w. 



(l-2Xi)Ja 



' v / (l + l/«)-(2A«-l) 2 



(Bll) 



Since, (1 + 1 /a) = (2A± — l) 2 and A_ < A,- < A + , the term 
under the square root in the above equation is always positive. 



1. The case: 8 < 1 

For any 8 < 1, there exists a C] given by the solution of 

1 — 2b + X+ = as 



e ci (8) = 



;i-5 2 r 



5 2 (3 + 10S + 35 2 )' 



(B12) 



and for 6 < 9 C1 the function g(X) has no singularities whereas 
it has one singularity for 9 > 9 C1 . As 8 — > we get 9 Ci ~ 
1 / (3 S 2 ) whereas 9 Cl ~ (1 - <5) 2 /4 as 8 -> 1". 

The 9 = 9 CI (8) line corresponds to the a = CC C] {8) line in 
the (a, 8) plane, where 



0c, (5) 



(l + 5)(l-5) 2 
3 + 10S + 35 2 



(B13) 
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2. The case: S > 1 

For 5 > 1, there again exists a 9 C2 given by the solution of 

1 - 2/?_A_ = as 



e C2 (S) = 



(8 2 -iy 

8 2 (3 + 105 + 35 2 )' 



(B14) 



and g(X) has either three or four singularities depending on 
whether > C2 or < Q C1 . In the limit 8 — s- °° we get 
9 C2 = 1/3 and 6 Cl — > as 8 — > 1. More precisely, 9 Cl ~ 1/3 - 
10/(95) as 5 -)• oo, whereas Q ~ (5- l) 2 /4 as 5 -> 1+ 

The = O (5) line corresponds to the a — a c - 2 (8) line in 
the (a, 8) plane, where 



a C2 (8) 



(l + 5)(5-l) 2 
3 + 105 + 35 2 



(B15) 



3. The case: 5 = 1 

It is instructive to illustrate the particular case of 5 = 1 , for 
which we have a — 9/2 and a = 0/4. Here from Eqs. ( fT7| i 
and d22l) we get 2b± = 9X± and 0A + A_ — - 1. It follows that: 



(a) 1 - 2/?+A_ = 2 and 1 - 2b+X+ = -9X+ < for > 0. 
This implies g(X) has a singularity at A — A a . We get 
A a = (l+2A+)/(2 + 0A+) andA a e (0,A+). 

(b) 1 - 2b_X+ = 2 and 1 - 2Z?_A_ = -0A_ > for > 0. 
This implies g(X) does not have any singularity at A = Aj,. 

(c) 1 + 2b+X + =2 + 0A+ > and 1 + 2b + X- = 0. However, 
since v(A_) =0, g(X) has a singularity at X = X c = A_. 

(d) 1 +2b_X+ = 2 + 0A_ > as 0A_ € (-1,0). Moreover, 
1 +2£>_A+ = and v(A + ) = 0. Therefore, g(X) has a 
singularity at A = Ad = A + . 

However, we have already seen that A* — > X± only when w — > 
=F°°. Therefore, for all practical purposes (any finite w) the 
singularities at X± are not relevant and hence we treat this case 
together with the case 5 < 1, > Cl where g(X) has only one 
singularity. However, for the 5 = 1 case, in principle, one can 
also use the results of Sec. |IV C| where the case of the three 
singularities is discussed. 



Appendix C: Steepest descent method with a branch point 



Let us consider the integral 



/ = 



2%i 



e rfuW 

V A a — X 



dX, 



(CI) 



where A a > 0. The position of the saddle point A* depends on 
the value of w, and depending on whether w > w* or w < w* 
we have A* < A a or A* > A a respectively. In the following, we 
consider the two cases one by one. 




I, 



Re X 



FIG. 9. (Color online) Schematic steepest descent contour (in red) 
for the case when the branch point is not between the origin and 
the saddle point X*. Here, it is shown for A* < 0, however, one 
can also have X* e (0, A a ). The direction towards which the contour 
bends, depends on the value of w. Here, it is shown for w < 0. For 
w > the contour bends towards the positive Re(A) axis, whereas 
for w = 0, the steepest descent contour is parallel to the Im(A) axis. 
The thick solid (orange) line along the Re(A) axis from A a represents 
the branch cut. 



1. The branch point is not between the origin and the saddle 
point: A a (0,A*) 

In this case, since A a lies outside the interval (0,A*), 
one can deform the contour of integration in Eq. ( |CT[ l into 
the steepest descent path through A* without hitting A a (see 
Fig. [9j. Along the steepest descent contour we define 

fw(X) - f w (X*) = -ic. 

Therefore, A a is mapped to a branch point at u 

and Eq. (jCTJ becomes 

qi (m) 



(C2) 
ib with 
(C3) 



with 



2ni 



qi(u) = g\{X) 



Vb- 



-. du 



y/b — iu dX 



VA a -A du 



(C4) 



(C5) 



Now, making a change of variable ^/ru — > u and taking the 
large-T limit we get 



2ni 



-?i(0)t 



-1/4 



where 



Vb dX 



: du, 



VAa-A* du 



(C6) 



(C7) 
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FIG. 10. (Color online) Schematic steepest descent contour for the 
case when the branch point A a is between the origin and the saddle 
point X*. The thick solid (orange) line along the Re(A) axis from A a 
represents the branch cut. The steepest descent contour goes around 
the branch cut as shown by C2 and C3 (in blue). The contribution 
coming from the circular contour (in magenta) C £ around the branch 
point becomes zero in the the limit of the radius e —¥ 0. The direction 
towards which the contours C\ and C4 (shown in red) bend, depends 
on the value of w. Here, it is shown for w > 0. For w < the C\ and 
C\ bend towards the negative Re(A) axis, whereas for w = 0, they 
are parallel to the Im(A) axis. 



Using -u 2 = i/"(A*)(A - A*) 2 + ■ ■ ■ as A -> A*, it can be 



found that 



dX 
du 



is/2 



(C8) 



Therefore, we get 



/ * oUA ) _^ ==Rl{Vxb)j (C9) 



where 



VX~^ s/2%xfl{X* 



: du. 



(CIO) 



We perform this integral in the Mathematica to get Eq. ( 30 1. 



2. The branch point is between the origin and the saddle 
point: A a e (0,1*) 

In this case, since A a lies in the interval (0,A*), the de- 
formed contour the through X* wraps around the branch cut. 
The contour of integration C = C\ + Co + C3 + C4 + C £ is 



shown in Fig. 10 The contour C £ represents the circular con- 



tour of radius e going around the branch point, and its con- 
tribution becomes zero in the limit e — > 0. The integral in 
Eq. (CI 1 can be written as / = +Ps, where Pq(w,t) is the 



contribution coming from the integrations along the contours 



C2 and C3, whereas P§ is the saddle point contribution com- 
ing from the integrations along the contours C\ and C4. In the 
following we evaluate Pg and P$. 



We consider, 



a. Branch cut contribution 



2ni Jc 2 +c 3 v A a - A 



dX. 



(Cll) 



We note that \JX A — A changes when one goes from C2 to C3. 
More precisely, A a — A = |A a — A|e'^, where (j) = +n on C2 
and (j) = —K on C3 (as (j) ~ for A < A a on the real-A line). 
Therefore, \/A a — A = +/|A a — A| 1//2 on C 2 and %/A a — A = 
— ; I Aa — A| 1,/2 on C3, using which from Eq. (Cll I we get 



Pb = 



1 



■ (A) 



|A-A a |V2 



dX. 



(C12) 



Since / W (A) is real, / w (A a ) > / W (A) > / W (A*) for A a < A < 
A*, and / W (A) is minimum at A* along the real A line, we set 



f w {X)-U{K) 



-2bu + u 



(C13) 



The branch point A a is mapped to u — 0. Using f' v (X*) = 0, 
we find that the saddle point is mapped to u = b, and b can 
be found by putting A = A* and u = b in the above equation, 
which gives Eq. (C3 I. With the above mapping from A to u, 



Eq. (CI 1 1 becomes 



ft = 



qi{u)- 



„—x{2bu—u 



- du, 



where 



«2(w) =£iW 



t/A 



From Eq. (C13 1, we get 

dA 
<t/m 



|A-A a |V2 du 
2{u-b) 



(C14) 



(C15) 



(C16) 



which is finite and nonzero everywhere between u = and 
u = b. Near u = we get 



c/A 



2/; 



-/;(Aa)' 



(C17) 



On the other hand, near u = b, by applying L' Hospital rule to 
Eq. ( |C16| l we get 



t/A 



V2 



(cis) 



Now, making a change of variable \/%u 
and then taking the large- T limit we get 



u in Eq. (C14i 



42 (0) 

T V4 



«2(Vt6), 



(C19) 
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where 



Hz) 



" \ e -2zu+u 2 

\fu 



du. 



(C20) 



The asymptotic forms of R 2 (z) can be easily determined from 
the above integral, which gives R 2 (z) ~ \f%j \/2~z as z^r oo. 
It can be shown that 



afA 



Therefore, 



| X — Aaj 1 / 2 a?M 



? 2 (0) = si (A,) 



c/A 



1/2 



(C21) 



dX 




du 


u=0 



1/2 



(C22) 



fa Saddle point contribution 



We evaluate this integral in Mathematica to get Eq. ( 37 1, 
where the generalized hypergeometric function has the series 
expansion 



2 F 2 (ai,a 2 ;b h b 2 ;z) = £ 



(ai) n {a 2 ) k Z n 



(C29) 



with (a) n = a(a + l)(a + 2) • • • (a + n — 1), (a)n = 1 being the 
the Pochhammer symbol. 

The large z behavior of R4 (z) can be found by expanding the 
term inside the square bracket in Eq. ( C28| > in powers of 1/z 
and integrating term by term. This gives ^4(2) ~ l/(2y/nz) 
for large z. 

On the other hand, # 4 (z) ~ Y(\/A)^/z/2n for small 
Z. Using this together with lim^*^ a y/b/\Xa — A*|'/ 2 
/"' (A * ) /2] 1 / 4 in Eq. \C2l\ we get 



We consider, 



1 r e T f»W 
Ps = ^-. gi(A) JA. 



(C23) 



We make a transform from A to m as defined by Eq. ( C2 1. In 
this case, the branch point A a is mapped to a branch point at 



// = ib where b is given by Eq. ( C3 1, and Eq. ( C23 1 becomes 



Ps = 



2ni 



q 3 (u) ■ 



with 



93(") =giW 



y/b + iu dX 
\JX a — X du 



du 



(C24) 



(C25) 



We found in the preceding sub-subsection that \/A a — A* = 
±;|A a — A*|'/ 2 below (+) and above (— ) the branch cut re- 
spectively. Therefore, <73(«) approaches two different limits 
as u — y form above (0 + ) and below (0~) respectively: 



93(0*) 



. gi(X*)Vb V2 
lAa-A*! 1 / 2 y/f^iX*)' 



(C26) 



where we have used Eq. ( |C8[ ) for the Jacobian. Thus, upon 
changing ^/tm — > u and taking the large- T limit yields 



ft 

where 

**(z) 



i (A* 



e~ M du f° e- 1 1 du 



(C27) 



y/z + i 



u r 

u J - 



- I due " 2 



Ti JO 



>/z + iu 
1 



\/z + iu \Jz — iu 



(C28) 



Ps 



r(l/4)gi(A*)e T ^ r ) 
2% [2r^(A*)]i/4 



(C30) 



Appendix D: Index of the notations 



7 is the viscous drag. 

k is the stiffness (spring constant) of the trap. 

D = kgT /y is the diffusion constant. 

Ty = y/k is the viscous relaxation time of the trap, which is 

introduced in Eq. ([TJ. 

To is the correlation time of trap modulation, which is intro- 
duced in Eq. Q. 
X c is defined in Eq. ( [TT] ). 
and 8 are given by Eq. 
a is defined by Eq. (B , and related to 6 and 5 by Eq. 
a is defined in Eq. ([^ 1. 



jU(A) is given by Eq. ( 1 1 
v(A) is given by Eq. ( 121 and Eq. (23 1. 
g(X) is given by Eq. (fl6|i. 



f w (X) is defined in Eq. (20 1. 



b± are defined by Eq. ( 17 1. 
A± are defined by Eq. (22 1. 
A is given by Eq. @. 

Xj with i £ {a,b,c,d}, are the positions of the branch points 
of g(X) and are given by Eqs. ( |B1[)- ( [B4| . 
h s (w) = f w (X*)is given by Eq. (25 1. 
ht(w) = f w {Xi) with i £ {a,b,c,dj~are given by Eqs. 
dE9V 



w*, with i £ {a,b,c,d}, are given by Eq. (B 
g(Xj) with i £ {a,b,c,d}, are defined by 
9 Cl , 9 C2 , and e c are given by Eqs. ( |B12) , ( |B14| i," and ((64] 
respectively. 



Hi}. 
Eq ([BIO 
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